Dynamical diffraction of high-energy electrons by light-atom structures: a multiple forward scattering interpretation

The T-matrix is used to compute the scattering of fast electrons by a regular array of effective spherical potential wells. An assessment of the forward scattering approximation and a real-space multiple scattering interpretation are provided.


Motivation
The theory of dynamical diffraction was first developed in 1928 (Bethe, 1928;Cowley, 1995), very shortly after the first experimental demonstration of electron diffraction. Over the 20th century, several theories of electron diffraction have been extensively applied in various areas of solid-state physics. In particular, the multiple scattering theory (MST) (Korringa, 1947(Korringa, , 1994Kohn & Rostoker, 1954;Dederichs, 1971) has proven to be very efficient for describing electronic properties of matter. Besides, MST provides an intuitive picture of dynamical diffraction. Partial wave scattering theory (Schiff, 1955) provides another description of multiple scattering with the possibility to take into account inelastic scattering (Howie, 1963) and loss of coherency (Howie, 2014). Other less commonly known attempts at an intuitive depiction of dynamical diffraction have been proposed to explain electron channelling in crystals, for example by VanDyck & Op de Beeck (1996). Dynamical diffraction must be taken into account when considering electron-atom interactions, since even at the very high electron energies commonly used in modern transmission electron microscopes, the interaction is so strong that, theoretically, the kinematic approximation is not valid for ideal crystals thicker than 10-15 nm (Hirsch et al., 1965;Glaeser & Downing, 1993;Subramanian et al., 2015). In practice, crystal growth cannot be controlled to such a degree of accuracy and, usually, nanocrystals of organic compounds have a size on the order of hundreds of nanometres. This is a challenging aspect of high-energy electron diffraction in crystallography, as it significantly complicates the structure determination process. However, structures have been successfully determined from ED (electron diffraction) data using the standard kinematic theory (Gemmi et al., 2019;Nannenga et al., 2014;Nannenga & Gonen, 2019). Although dynamical refinement (Palatinus et al., 2013) usually leads to better intensity predictions (Gemmi et al., 2019), the agreement between theory and experiment is still significantly worse than that obtained for X-ray data (Oleynikov et al., 2007;Klar et al., 2021).
Although there is no quantitative prediction of the effects of dynamical diffraction on the ability to solve a crystal structure, it may be speculated that dynamical diffraction may be attenuated by other effects such as a lack of coherency caused by solvent scattering or crystal defects. Therefore, more accurate theoretical models are needed for getting better results from electron diffraction experiments.

State of the art
The theory of multiple scattering has a long and rich history, including modern theoretical developments and implementations (Sé billeau et al., 2017;Zabloudil et al., 2005). However, in the context of high-energy electron diffraction, the multislice (MS) (Cowley & Moodie, 1957) approach and Bloch-wave (BW) (Bethe, 1928;Metherell & Fisher, 1969) approach are the most popular methods for simulating diffraction patterns by crystals.
MS is particularly well suited for solving large problems as it involves successive convolutions, which can be very efficiently computed with the fast Fourier transform (FFT) (Ishizuka & Uyeda, 1977). However, in order to avoid aliasing, transverse periodic boundary conditions must be met. This is not always possible if the crystal is in an arbitrary orientation as is the case in typical continuous-rotation electron diffraction experiments. A small beam tilt can also be used as an option to remedy this aspect but only extends reliably to a maximum 3-6 (Ishizuka, 1982;Chen et al., 1997). Simulations with arbitrary orientations could also be performed by considering all unit cells in a crystal, while adding zero padding at its ends. Even though modelling a full crystal can quickly become intractable computationally, the use of a smooth envelope function may improve this commonly admitted limitation of MS (Kirkland, 2019). A few very efficient MS implementations have been reported to date (Ophus, 2017). Some of them are targeted at convergent-beam electron diffraction (CBED), while others have special features such as inclusion of inelastic scattering (Allen et al., 2015). Yet, an approach combining full modelling capabilities with an acceptable computational efficiency, designed specifically for the continuous-rotation electron diffraction of large organic structures, does not seem to be available. Such an approach would indeed be of significant value for macromolecular structure determination from ED data.
On the other hand, the BW method can simulate ED for structures in arbitrary orientations but does not scale well with the structure size. Although a few rather efficient BW implementations (Zuo & Weickenmeier, 1995) have been proposed to include a large number of beams, BW can hardly be applied to very large structures since the scattering matrix becomes prohibitively large. Due to the unfavourable OðN 3 Þ complexity of matrix diagonalization, where N is the number of beams, the method is hardly applicable in practice to the determination of macromolecular structures where it would have to be done numerous times. Non-periodic structures, defects and solvent scattering are also hard to model with this method.
In other areas of physics, specialized multiple scattering approaches have also been developed. In particular, the T-matrix approach has been extensively applied to wave propagation related problems such as electromagnetics (Hamid et al., 1990a,b;Eremin et al., 1995), optics (Moine & Stout, 2005) and acoustics (Silva et al., 2012;Godin, 2011). Although it does not compete with MS from a computational standpoint, this approach benefits from providing an exact solution to Schrö dinger's equation for an ensemble of spherically symmetric atomic potentials in the independent atom model (IAM) approximation. Besides, ab initio realspace multiple scattering calculations have also been developed in X-ray photon emission spectroscopy with elaborate potentials including many-body effects (Rehr et al., 2009). Similar approaches based on the T-matrix also exist in electron energy-loss spectroscopy (Sé billeau et al., 2006).

Contribution and outline
The purpose of this paper is to adapt the T-matrix formalism specifically to the case of scattering of fast electrons by lightatom structures. In particular, we present an intuitive picture of multiple scattering in the forward scattering approximation. We propose a comparison with the MS interpretation of multiple scattering. In addition, we discuss the validity of the forward scattering approximation and the phase grating approximation used in the MS approach. Finally, we draw conclusions and outline extensions of the proposed approach in order to account for incoherent inelastic electron scattering.

Theory
The scattering of fast electrons by an atomic structure is described by the wavefunction É, obeying the following Schrö dinger's equation (Kirkland, 2019): where his Planck's constant, m e the mass of the electron, e the elementary charge, VðrÞ the spatially varying electrostatic potential created by the atoms, and E the electron energy.
In an ED experiment, E ) Vðr ! 1Þ ¼ 0, which corresponds to the continuum, non-quantized state of the system (Chuang, 1996). Therefore, we solve (1) with E as a parameter, equal to the energy of incident electrons.
The total wavefunction É is described as the sum of an incident wave É ðiÞ and a scattered wave É ðsÞ : The incident wave is assumed to be known, i.e. a plane wave or a Gaussian wave for example. The T-matrix aims at determining the scattered wave É ðsÞ by the sample.

T-matrix formulation
In its standard form, the T-matrix approach solves for the case where the incident wave is described by a plane wave of wavenumber k 0 ¼ 2= (optics convention), and the electrostatic potential is modelled by a uniform constant inside nonoverlapping spheres. The constant potential and radius of the spheres depend on the atom type.
Although this is not an accurate representation of the actual potential profile VðrÞ, using an appropriate amplitude for the constant potential V 0 should provide similar values of scattering cross sections. The scattering cross section being a reliable figure of merit of the strength of an interaction, one can assume that a suitable choice of radius and constant potential should therefore provide a faithful depiction of the extent of multiple scattering in the structures of interest.
The setup is shown in Fig. 1. This formulation is well established (Hamid et al., 1990a,b;Eremin et al., 1995) and the theory is outlined here for the purpose of introducing the forward multiple scattering approximation. In the T-matrix approach, the electrostatic potential is assumed constant inside different regions of the 3D space. In each region, the problem is reduced to the Helmholtz equation in spherical coordinates: where V p ! 0 is the constant positive potential inside sphere p of radius a p centred at d p , k p the wavenumber inside the sphere and n p is sometimes referred to as the refractive index by analogy with optics. In the remainder, n p will be referred to as the potential strength for clarity purposes. The scattered wavefunction É ðsÞ can be decomposed into the scattered wavefunctions of individual atom spheres, where N is the number of spheres.
where p ¼ 1 . . . N, j l and h ð1Þ l are the spherical Bessel and Hankel functions of the first kind, respectively, Y m l are the spherical harmonics of order l and azimuthal order m. Note that these equations are expressed in the reference frame of each sphere p, hence the use of variable r p .
The unknown coefficients a p;lm , b p;lm are found by imposing the continuity of the wavefunction and its radial derivative at the surface of each sphere p. After some mathematical manipulations (see Appendix A), this results in the following linear system of equations: where I is the identity matrix, A the unknown vector of coefficients, T the cross-coupling matrix and L the known incident wave which appears on the right-hand side of (6). Note that, in this formulation, the unknown coefficients, found from the continuity interface conditions, will model the response of each individual sphere to both the incident wave and the waves scattered by the other atoms. This model therefore does account for multiple scattering.

Far field and scattering cross section
In electron crystallography, diffraction patterns are recorded that correspond to the intensity of the scattering amplitude profile jf ð; Þj 2 for various crystal orientations. A diffraction pattern is obtained in the far-field diffraction 182 Tarik   Setup and conventions for the scattering problem solved with the Tmatrix approach. A plane wave with a wavevector k 0 with incident angles ( i , i ) interacts with N spherical atoms, each with individual location (for example, atom p at d p ), radius (a p ) and electrostatic potential strength (n p ). The scattered wave is observed at observation point P at r in the global reference frame and at individual distances in the local reference frame of the spheres (for example, jr p j). See discussion in the text for additional details. In practice, the incident plane wave is taken along the +z axis and the structure oriented accordingly. regime, which arises when the shape of the angular radiation scattering amplitude f ð; Þ no longer depends on the distance to the scattering specimen. Using the asymptotic behaviour h ð1Þ l ðk 0 r p Þ ' ðÀjÞ lþ1 f½expðjk 0 r p Þ=ðk 0 r p Þg and since p ' , p ' , the far-field scattering amplitude f p ð; Þ from each individual sphere p can be written as The total far-field scattering amplitude f ð; Þ is the sum of the scattered field contributions from all individual spheres. Note that the scattered field from each individual sphere already accounts for the scattering between the spheres through the coefficients a p and b p . Since, in the far field, The normalized differential scattering cross section can be obtained as where we have used É ðsÞ ðr; ; Þ ' r!1 ½expðjk 0 rÞ=ðk 0 rÞf ð; Þ.
3. Real-space forward multiple scattering picture

T-matrix forward scattering approximations
Equation (6) is a convenient way to represent the system as it readily identifies L as the solution to the uncoupled (T = 0) problem. If T = 0, the system can be broken down separately into the scattering problems for each individual sphere, for which the solutions are immediately identified as a p;lm ¼ c lm u p;l , b p;lm ¼ c lm v p;l . These are indeed the well known analytical solutions of Mie scattering by a soft sphere (Balanis, 1989).
The cross-coupling matrix T accounts for multiple scattering effects. If A is written as A ¼ ð. . . a p;lm ; b p;lm . . .Þ T , then T is a matrix with only off-diagonal components: where T pq represents the scattering from sphere p due to the scattering from sphere q. By considering all components T pq , T qp , we assume that the scattering from sphere p affects scattering from sphere q and vice versa. Matrix (10) is completely filled with nonzero entries apart from its diagonal components. It is therefore not a sparse matrix, which implies quadratically growing memory requirements. Moreover, since full intersphere scattering is considered, the solution of equation (6) requires matrix inversion which is computationally expensive. However, as detailed below in Section 4.2, backscattering can be neglected for very fast electrons, which is known as the forward scattering approximation. This results in T being lower triangular if the spheres are sorted in ascending order along e z . In the case of atoms lying in the same coordinate plane, i.e. slice, the atoms can be sorted according to their transverse coordinates. Then, as 90 scattering is neglected in the forward scattering approximation, the T-matrix remains lower triangular. Inversion of the lower triangular matrix is computationally tractable and calculations can be performed sequentially one slice after the other. Therefore, the forward scattering approximation converts the implicit self-consistent T-matrix scheme into a forward scattering explicit scheme similar to the approach taken in the MS method.

Multiple scattering approximations
Since A 0 ¼ L represents single scattering, we can establish that A 1 ¼ TA 0 accounts for secondary scattering. Similarly, outward scattering amplitudes from sphere p can be written as where the first term accounts for kinematic scattering, the second term for secondary scattering, the third term for tertiary scattering, and so on. This is a development similar to the Korringa-Kohn-Rostoker (KKR) theory of multiple scattering (Korringa, 1947(Korringa, , 1994Kohn & Rostoker, 1954). This forward multiple scattering picture is illustrated in Fig. 2 for the case of a two-atom system. From a computational point of view, equation (11)  Two-level forward scattering approximation using the T-matrix approach with an incident plane wave propagating upwards. The third column represents scattering from the top scatterer, located at distance kd from the bottom scatterer. By neglecting backscattering, this can be approximated by the sum of (i) the kinematic scattering of the incident plane wave by the top scatterer (first column) and (ii) the secondary scattering following kinematic scattering from the bottom scatterer (second column). The scattering from the bottom scatterer is viewed by the top scatterer as a particular incident wave (neither planar nor spherical) emitted from a point source at a distance kd. expansion depends on precomputed coefficients b ð0Þ p;lm , computation of the b ðnÞ p;lm coefficients for the nth scattering term can be massively parallelized for each atom.

Scattering probabilities
The probability of an electron being scattered elastically n times can be estimated classically from ballistic arguments using a continuous model of matter. It can be established (Egerton, 2011) that the probability of an electron being elastically scattered n times after passing through an amorphous sample of thickness z follows a Poisson distribution P n e ¼ ð1=n!Þðz=l e Þ n expðÀz=l e Þ, where l e is the elastic mean free path. It is paramount to note that this classical expression models incoherent scattering, not including interference effects between multiple scattering events. This model is therefore valid for classical particles not exhibiting particlewave duality behaviour. In quantum systems, it is also possible to consider incoherent scattering if inelastic scattering is strongly present. Indeed, inelastic scattering not only induces a transfer of energy between scatterers but also a loss of coherency through phase shifts (Howie, 2014). As a result, quantum systems, where inelastic scattering is predominant, would scatter incoherently. For example, inelastic scattering may be due to electron-electron interactions, solvent or impurity scattering, as well as thermal diffuse scattering.
In this section, we suggest a way to consider probabilities of multiple coherent elastic scattering events. It will be shown that while multiple scattering amplitudes can be considered, it may not be possible to identify corresponding probabilities without any ambiguity precisely because of interference.
Using a wave-like description, the probability of an electron undergoing a scattering event can be determined from the scattering cross section as follows. Let S be the area over which the specimen is illuminated.
Let z be the thickness of the specimen. The incident electron plane wave can be normalized by jAj 2 ¼ 1=Sz so that it integrates over the interaction volume V ¼ Sz to a single electron. Then, the flow of electrons per unit time per unit area (Vainshtein, 1964) By definition of the cross section, J 0 is the probability of an electron being scattered every second. Since it takes a time Át ¼ z=v 0 for an electron to pass through the specimen, the overall probability of a single electron being scattered is therefore P scat ¼ J 0 Át. Therefore, after inserting the expressions for J 0 and Át, where S can for example be the area of the aperture in selected-area electron diffraction (SAED). When performing a MS simulation, S is the area of the simulated domain (transverse supercell if performing a periodic simulation).
Here, we define f n as the complex scattering amplitudes obtained by considering scattering from the wave being scattered n times. It can readily be established that We can then define P n ¼ n =S as the probability of an electron being scattered n times. The choice of n-time scattering cross section n is arbitrary precisely because of the cross-coupling interference terms f n f Ã m . It has been chosen in such a way that when more n-scattering terms are taken into account, probabilities of higher n-level scattering events are transferred from the lower n-level scattering events. In the case of an application to two-beam scattering theory, the intensity of the forward and scattered beams oscillates as the thickness of the sample grows. The scattered beam is zero when the forward beam has been scattered an even number of times and maximum when the forward beam has been scattered an odd number of times (Cowley & Moodie, 1957). Our definition of scattering probabilities should reflect this behaviour.
The probability of an electron not being scattered is naturally P coh ¼ 1 À P scat , where P scat ¼ P n P n . Note that, since S ) , the probability of scattering is necessarily less than 1.
For a very large number of scatterers, N, regularly spaced by distance dz, we can define the average cross section a ¼ =N.
The corresponding scattering probability is then PðdzÞ = a dz ¼ dz=l e where ¼ 1=Sdz is the density (since Sdz contains only one scatterer) and l e ¼ 1= a is the mean free path. This is consistent with the definition used by Latychevskaia & Abrahams (2019).

First Born approximation and kinematic scattering
Let us assume normal incidence with i ¼ 0, so that e p ¼ expðjk 0 d p Þ, and note that the Ewald sphere can be represented by q = Kðsin cos ; sin sin ; 1 À cos Þ in reciprocal space (or transfer vector momentum space), with K ¼ 1= being the wavenumber in crystallographic convention. Keeping only the kinematic term in equation (11) and inserting in the far field in equation (8) give where f ðeÞ p is the scattering amplitude of sphere p in the first Born approximation, also known as the form factor. Equation (13) corresponds to the standard kinematic formula for the structure factor traditionally used in crystallography. Note that the kinematic approximation is nothing else than the first Born approximation applied to the whole assembly of atoms.
It will be shown in Section 4.2 that the first Born approximation holds in the case of single scatterers, but fails in application to multi-atom systems. This means that the kinematic approximation is not generally valid.

Multiple scattering in the MS approach
In MS, the forward scattering approximation is used and the potential discretized in slices. The wavefunction is propagated from one slice to the other using a Fresnel propagator.
A multiple scattering interpretation of the MS approach has been proposed (Cowley & Moodie, 1957), which, although analogous to the one presented above, differs in that it is stated in reciprocal space. More details on the differences between the two interpretations can be found in Appendix B.

Application and results
Very efficient open-source packages exist for T-matrix calculations, but they are targeted at applications in electromagnetics (Egel et al., 2017;Kottke, 2020). There are also packages targeted at photon and electron spectroscopy and related to the T-matrix, such as msSpec (Sé billeau et al., 2011) and FEFF9 (Rehr et al., 2010).
Here, an open-source package, adapted to high-energy electron diffraction and including the forward scattering approximation, has been developed and made available (Drevon, 2021).

Validity of the implementation
In practice, the size of the matrix in equation (6) has to be truncated to a maximum integer order max due to computer memory limitations.
A rule for obtaining accurate results is to use an integer larger than the maximum value of the normalized radius ka, which we call ka max ¼ max p ka p where a p is the radius of the spheres and k the wavenumber.
Besides, the translational addition theorem (Dufva et al., 2008) used for expressing the scattered field of any sphere in the reference frame of another sphere introduces an approximation of the translated spherical Hankel functions. This theorem is also more commonly known as the structure constant expansion or Kasterin expansion. The computational accuracy of this translation decreases from the centre at which this expansion is written. This is analogous to a Taylor expansion, for which the accuracy away from the point at which it is expressed increases with increasing expansion order. Note that the distances between the spheres do not affect the accuracy of this translation operation and should not affect the choice of max . Once again, the radius of the spheres determines the choice of max . This is illustrated in Fig. 3(a), where the error between h ð1Þ l Y 42 , computed at the origin, using the translational addition theorem with d p ¼ ð0; 3; 5Þ and max ¼ 10, is displayed in log scale. The choice of max is best evaluated by assessing the continuity of the wavefunction at the surface of the spheres as shown in Fig. 3(b) for N ¼ 4 and ka max ¼ 4. This is also used to validate the correctness of the implementation since it can be seen that machine accuracy can be reached when increasing the order of the expansion.

Validity of forward scattering and phase grating approximations for light atoms
In the case of very fast electrons typically used in transmission electron microscopes, E = 50-300 keV. Inclusion of relativistic effects results in a wavelength = 0.02508 Å at 200 keV, which will be assumed from now on unless stated otherwise.
In the IAM approximation, the Coulomb potential is created by the charge of the nucleus and its electron cloud. It is typically fitted with a sum of three screened Coulomb potentials and three Gaussian terms (Kirkland, 2019) which results in the potential shown in Fig. 4(a).
The solution to Schrö dinger's equation in such a potential can only be solved perturbatively (Mü ller, 1965) and is beyond the scope of this study. However, indicative values for k p and ka p can be used with a multi-shell representation as shown in the form of blue patches in Fig. 4(a). Although the range of the     screened Coulomb potential is theoretically infinite, it can be truncated to radius ka for most practical purposes. This Coulomb potential truncation is commonly known as the muffin-tin model. The number of concentric spheres to use has mainly an effect on the large-angle representation of the form factor which does not impose severe constraints in a low-angle forward scattering approximation. In Fig. 4(b), the multi-shell scattering amplitudes calculated in the Born approximation are shown for increasing values of truncation radius. The figure only shows a satisfactory agreement with the electron diffraction scattering factors for normalized radius as large as ka ¼ 400. Even though the potential VðrÞ is extremely small at such a large value of the radius ka ¼ 400, those shells are required to account for the proper low-angle representation of the form factor. As mentioned above, ED simulations based on the T-matrix approach are very expensive computationally at ka ¼ 400, because higher-order terms of series in equations (4), (5) would need to be included. However, the medium truncation range should be sufficient for providing a reasonable picture of dynamical scattering. Fig. 5(a) shows the total scattering cross section of a single sphere with increasing radius for a range of values of potential strength n p . The dark blue curve shows the locations of the spherical shell for carbon. It becomes almost flat for the parameter set (ka ¼ 30, n p ¼ 1:001), from which it keeps increasing very slightly to reach the asymptotic value ' 0:0035. This value is almost identical to the average elastic cross section for real carbon in the Born approximation (Latychevskaia & Abrahams, 2019) which confirms that the chosen parameters capture the strength of incident electron-carbon interaction. Therefore, even though the chosen parameter set is a coarse representation of the potential profile of a real carbon atom, it should provide a reasonably faithful estimate of the extent of multiple scattering of real carbon-based samples, which is the main objective of this study. Fig. 5(b) shows the shape of the far-field amplitudes for two sets of ka, n p values. The first one is for the extreme value of ka ¼ 3 and n p ¼ 1:1 while the second corresponds to one of the carbon atoms. While differences are obvious for the first case, the Born approximation is very accurate for the carbon atoms, although some minor differences appear at low angles where the phase grating approximation provides some improvements too. Fig. 6(a) shows a ðka; n p Þ map of the difference between the exact coefficients of the two-body problem and those calculated using the kinematic approximation. This error is calculated as errðb p;l;m Þ ¼ P p;l;m jb p;l;m À b approx p;l;m j: ð14Þ Fig. 6(b) shows the same thing as a result of using the forward approximations.
For this case, a value of kd ¼ 3ka is used, which is close to interatomic distances. Overall, it is clear that the forward scattering approximation is very accurate over all ranges of parameter sets for carbon atoms, which is expected since there is very little backscattering beyond 90 . On the other hand, the kinematic approximation does not appear quite as good, even for this, a mere two-atom problem. Although not shown here, both the kinematic and forward scattering approximations tend to work slightly better with increasing distances kd, since coupling between the spheres gets reduced and, therefore, is less likely to affect scattering from the other spheres. Using low values of n p results in an overall good approximation of both the uncoupled and forward scattering approximation. This is an anticipated result, since for weak potentials the kinematic approximation works better. The uncoupled approximation improves with decreasing radii, since a small ka results in a small scattering cross section. On the other hand, the forward scattering approximation improves with larger values of ka, which make backward scattering less likely. Fig. 7(a) shows the scattering coefficients b p;l;m in complex space for the given two-body problem, where n p ¼ 1:01, kd ¼ 3ka, while increasing radius ka (up to ka ¼ 15) for a 186 Tarik    Error of b p;lm [equation (14)] in ðka; n p Þ space (colour axis in log scale) for a two-scatterer system using the (a) kinematic scattering approximation, (b) forward scattering approximation. The blue dots correspond to the location of the spherical shells of a carbon atom at E = 200 keV (shown only for first five shells for visualization purposes). selected set of orders l. Only the coefficients for sphere p ¼ 2 are shown since it is the sphere most affected by the approximations. Here, only m ¼ 0 coefficients are used, because in the case of planar illumination, this configuration has azimuthal symmetry. The coefficients are computed using the full Tmatrix, the forward scattering approximation and the kinematic approximation. In order to get this picture, the coefficients were divided by the phase factor at sphere p, i.e. b p;lm expðjkdÞ. The arms of the spiral are rotated by =2 from one order to the other as a result of the j l factor in the spherical expansion of the incident plane wave. All the coefficients increase with increasing ka as a result of the increasing scattering cross section with increasing ka. For a given ka, the strongest contributing term orders are around l ¼ ka as a direct manifestation of the convergence behaviour presented above. It can be observed that the kinematic approximation error increases with the radius while the forward approximation is very accurate across all radii. Fig. 7(b) shows a comparison of the computed far-field scattering amplitudes for N ¼ 2, ka ¼ 30, n p ¼ 1:01, kd ¼ 2ka using the full T-matrix, the forward approximation and the kinematic approximation. The peaks and valleys are both due to the single scattering profile of the constant sphere and the interferometric path length expðjkdÞ. It is possible to see that the coupled problem has a slight averaging effect over the kinematic pattern, a well known feature of dynamical diffraction.

Successive multiple scattering approximations
Finally, we consider an array of N identical scatterers regularly spaced by kd ¼ 3ka which correlates with lengths in molecules.
Here, we consider the successive multiple scattering approach under normal illumination. This is an extreme case where strong dynamical scattering is expected. (14)] on scattering amplitude coefficients for an increasing number of spheres while also varying the number of successive approximations while keeping parameters ka ¼ 7, kdka ¼ 3 fixed.

Figs. 8(a) and 8(b) show the evolution of the error [equation
The forward scattering approximation becomes worse for increasing numbers of spheres as a result of accumulated errors. This is a limitation of a forward propagation scheme such as the one used in the MS approach, whose accuracy gets worse with increasing simulated sample thickness (Kirkland, 2019). Overall, the successive approximation schemes naturally converge to the forward approximation for sufficiently large n. However, higher-order terms need to be included in equations (4), (5) as the number of spheres (or thickness of the sample) increases. This result is similar to other multiple scattering approaches. For example, in the case of a two-beam configuration, it was shown (Cowley & Moodie, 1957) that the scattered beam could be expressed using expansion (17), considering only the scattering terms of the primary beam n ¼ 1; 3; 5 . . .. The famous two-beam analytical scattering expression is obtained as the sum of the infinite series. Fig. 9(a) shows the scattering probabilities of n-times scattering for ka ¼ 7, n p ¼ 1:01, ka ¼ 3kd with increasing N. The dynamic scattering probability P dyn ¼ P n>1 P n and the kinematic as well as the total scattering probabilities are shown. It seems that all multiple scattering terms increase with the number of spheres, in contrast to the ballistic picture of the Poisson distribution (Egerton, 2011), also shown as dasheddotted lines calculated using the average scattering cross section. The main reason for these trends is that the ballistic picture ignores the effect of interference, which drastically affects the overall scattering process. Fig. 9(b) shows the scattering amplitudes f n , associated with n-times scattering for ka ¼ 7, n p ¼ 1:01 and N ¼ 15. The scattering amplitudes have similar shapes, which mostly contribute to constructive interference, although slight peak misalignments can be seen. The reason for the large-angle   (14)] on the scattering amplitude coefficients with increasing number of spheres N for ka ¼ 7, n p ¼ 1:001, kdka ¼ 3. Squares: forward scattering approximation; diamonds: kinematic approximation; triangles: secondary approximations; coloured dashed lines: successive n-times multiple scattering approximations. (b) Error on the scattering amplitude coefficients in the increasing successive approximations for ka ¼ 7, kdka ¼ 3 and N ¼ 100 spheres while varying n p .

Figure 7
(a) Complex space scattering coefficients b p;l;m for a two-body problem for n p ¼ 1:01, kd ¼ 3ka with increasing radius ka. Only the coefficients for sphere p ¼ 2 and m ¼ 0 are shown. The coefficients were computed from full T-matrix (open circles), in the forward scattering approximation (squares) and the kinematic approximation (diamonds). (b) Comparison of far-field scattering amplitudes for N ¼ 2, ka ¼ 30, n p ¼ 1:01, kd ¼ 2ka using full T-matrix (solid blue), in forward approximation (green dashed line) and kinematic approximation (red dashed-dotted). The forward approximation is indistinguishable from the exact solution in this case.
scattering, visible on the figure, is due to the small radius used in this example. As mentioned above, a much larger normalized radius would be required to properly account for the lowangle representation of the far-field diffraction pattern. However, because of the comparable scattering cross section of the chosen model and a real carbon atom, the extent of multiple scattering should have been reasonably estimated.
Note that the examples presented here have been performed for electron incident energies of 200 keV, which is a typical value in transmission electron microscopy. Decreasing the energy would result in lower normalized radii, which would correspond to shifting the blue curve for carbon in Fig.  5(a) to the left. As a result, spherical harmonics of lower orders would be sufficient to achieve good accuracy. On the other hand, the forward scattering approximation would lose accuracy at smaller thicknesses, and then full inversion of the system (6) might be required. For higher energies, the forward approximation becomes more accurate but the size of the system increases, making the problem computationally harder.

Conclusion and perspective
An alternative approach, based on the T-matrix formalism, has been applied to the scattering of fast electrons, by lightatom structures. The validity of important approximations, used in the MS approach, has been discussed, and a multiple scattering approximation framework has been proposed and put in perspective compared with other existing interpretations. While it is worth mentioning that the T-matrix does not scale favourably with large increasing values of the wavelength-normalized radii ka, as commonly used for high-energy electron diffraction, it can still provide valuable insights in understanding dynamical diffraction effects thanks to its ability to solve the wave equation exactly. For a 200 keV electron beam, each spherical coefficient a p , b p accounts for one wavelength of physical space comparable with the sampling resolution used in typical MS simulations. For a truncation longitudinal number n ¼ 400, the total number of coefficients for all atoms would be 2n 2 N since azimuthal coefficients m 6 ¼ 0 are necessary for nonlinear array arrangements. For N = 10 000 atoms, this is on the order of a few gigabytes of memory. While inverting a matrix of this size would be computationally intractable, the forward scheme approximation, presented in this paper, makes such a computation viable. Indeed, it offers the possibility of parallel computation for all atoms in the same slice, similar to the way it is traditionally done in finite difference schemes.
In the presented approach, the spherically symmetric effective potential does not model the real electrostatic potential very accurately. However, we anticipate that the proposed interpretation of multiple scattering is equally applicable to the more accurate case of a screened Coulomb potential. Such a potential could naturally be included by using a basis of radial functions representing solutions to the homogeneous Schrö dinger's equation, which can be obtained numerically. In fact, since the T-matrix approach mainly relies on the spherical symmetry of the individual scatterers, it can be adapted to any family of radial functions, provided that translational coefficients (15), (16) can be computed numerically.
For periodic structures, the Green's functions (GFs) could be used as basis functions. All GFs are solutions to Schrodinger's equation in a periodic potential, expressed as an expansion upon the lattice vectors of the crystal under consideration. They would automatically satisfy the boundary conditions imposed by the Bloch theorem in crystals. This is the starting point of the Kohn-Rostoker method (Kohn & Rostoker, 1954).
An apparent limitation of both T-matrix and MS approaches lies in the use of the IAM, which, by definition, ignores the effect of bonding. Theoretically, such bonding could be included (Natoli et al., 1986). In fact, the MS approach can be adapted, at a computational cost, to bonding models such as the transferable aspherical atom model (TAAM) or even potentials determined by density functional theory. On the other hand, the single-site T-matrix approach, presented in this paper, cannot work directly with such potentials because the spheres are treated as non-overlapping. A possible option would be to employ effective potentials calculated in the framework of self-consistent electronic structures, where the effect of bonding manifests itself in the bonding-corrected effective potential. However, it is still an open question whether the effects of such bonding play an important role in structure determination of organic structures by electron diffraction.
The advantage of the multiple scattering approximation is that it may offer the possibility to include consistently the effect of incoherent and inelastic scattering, which can provide a greater insight into dynamical diffraction effects in real experiments. Incoherent scattering could, for example, be included in a stochastic fashion as randomly affecting the phase of the scattered waves. Inelastic scattering might require the use of dedicated electrostatic absorption potentials. It is 188 Tarik  (a) Scattering probabilities (dashed coloured lines) of n-times scattering, kinematic (cyan solid line), dynamic scattering (blue solid line) and total scattering (black solid) probabilities. The total scattering probability following the Poisson distribution using the average scattering cross section is shown as the black dashed-dotted line. Probabilities have been integrated over angular distribution using ka ¼ 7, n p ¼ 1:01 and ka ¼ 3kd. (b) Scattering amplitudes f n , associated with n-times scattering (for n = 1-, 2-and 3-level scattering) for ka ¼ 7, n p ¼ 1:01 and N ¼ 15. The total scattering amplitude f ¼ P n f n is shown in solid black.
anticipated that incoherent and inelastic scattering may have a significant effect in dynamical diffraction even when energy filters are used.

Data availability statement
Computer codes, developed in this study, are freely available under the GNU GPL licence v3.0 from https://github.com/ ronandrevon/pyScatSpheres. No new ED data were produced in this study.

APPENDIX A T-matrix
The unknown coefficients a p;lm , b p;lm are found by imposing the continuity of the wavefunction and its radial derivative at the surface of each sphere p: where É ðiÞ is the incident electron wavefunction, É ðinÞ p the scattered wavefunction of sphere p expressed inside the sphere, É ðoutÞ p the scattered wavefunction of sphere p expressed outside sphere p, and a p the radius of sphere p.
Using the orthogonality of the spherical harmonics, the following linear system yields the unknown coefficients: a p;lm ¼ u p;l c lm þ u p;l P N where the translational addition theorem (Dufva et al., 2008) has been used to express the field, scattered by sphere q in the reference frame of sphere p, formally written as É ðoutÞ q ðr p Þ. This operation involves the coupling coefficients a ðout-inÞ ;;l;m ðd pq Þ, where d pq ¼ d q À d p .
The coefficients c lm are related to the incident wave. In the case of a plane wave expðjk 0 Á rÞ, j ¼ ffiffiffiffiffiffi À1 p , the addition theorem is used to expand the plane wave on the spherical Bessel solutions basis set: Þe p e p ¼ expðjk 0 d p p Þ p ¼ sinðÂ p Þ sinðÈ p Þ sinð i Þ þ cosðÂ p Þ cosð i Þ; where d p ; Â p ; È p are the spherical coordinates of the centre of sphere p in the global coordinate system and 0 i is the angle of incidence with respect to the e z axis. For propagation in the ðy; zÞ plane, assumed in this paper, i ¼ =2 and e p is the phase offset at sphere p. The different notations are illustrated in Fig. 1. The coefficients u p;l and v p;l are expressed as u p;l ¼ h 0 l ðk 0 a p Þj l ðk 0 a p Þ À h l ðk 0 a p Þj l ðk 0 a p Þ 0 j l ðk p a p Þh 0 l ðk 0 a p Þ À n p j 0 l ðk p ad p Þh l ðk 0 a p Þ v p;l ¼ n p j 0 l ðk p a p Þj l ðk 0 a p Þ À j l ðk p a p Þj 0 l ðk 0 a p Þ j l ðk p a p Þh 0 l ðk 0 a p Þ À n p j 0 l ðk p a p Þh l ðk 0 a p Þ ; where z 0 l ¼ @ z l ðÞ. Each of equations (15), (16) represent n 2 equations for individual spheres where n ¼ l max is the truncation longitudinal order. Note that each equation involves the outside coefficients b q;l;m of all the other spheres through the translational coefficients. As a result, equations (15), (16) represent a N Â 2 Â n 2 system with as many unknowns as the number of equations and, therefore, can be written in the matrix form (6).

APPENDIX B Multiple scattering picture of the MS approach
The expression for the scattering amplitude f ðh; kÞ of beam h; k is proportional to the following expression: where h ¼ P N n¼1 h n , k ¼ P N n¼1 k n , h u ¼ P NÀ1 n¼1 h n , k v ¼ P NÀ1 n¼1 k n , l w ¼ P NÀ1 n¼1 l n , H ¼ NÁz is the total thickness of the sample, N the number of slices of thickness Áz, Q h;k;l = Àj=Áz h;k expðÀ2jl n z n =cÞ + F h;k;l , F h;k;l is the structure factor, ¼ 2m e h= the interaction parameter, the excitation error of beam ðh; k; lÞ and n is the excitation error of beam ð P n r¼1 h r ; P n r¼1 l r Þ. The excitation error is expressed as where a, b and c are the lattice constants of the crystal and K ¼ 1= the wavenumber. Equation (18) expresses the longitudinal distance of beam ðh; k; lÞ to the paraboloid shown in Fig. 10(a). This paraboloid is a very accurate representation of the Ewald sphere for large wavenumber K. Therefore, is very close to the excitation error commonly defined in crystallography. From equation (17) it is possible to gather terms in powers of Q NÀn 0;0;0 into f ðnÞ h;k , which corresponds to the incident beam scattered n times before contributing to reflections h; k. The term Q N 0;0;0 only appears for h ¼ k ¼ 0 and corresponds to the unscattered direct beam. There are N terms involving the factor Q NÀ1 0;0;0 depending on which slice the single scattering research papers Acta Cryst. (2023). A79, 180-191 event took place. There are NðN À 1Þ=2 terms involving the factor Q NÀ2 0;0;0 depending on which pair of slices the two-level scattering process is considered, and so on for the higher-order multiple scattering terms (Cowley & Moodie, 1957). This multiple scattering picture is illustrated in 2D in Fig. 10(b) for a case with N ¼ 3 and using only beams in the zero-order Laue zone (ZOLZ) l i ¼ 0.
Using only ZOLZ beams, f